Stochastic Turing patterns in the Brusselator model 



Tommaso Biancalani 
Dipartimento di Fisica, Universita degli Studi di Firenze, 
via G. Sansone 1, 50019 Sesto Fiorentino, Florence, Italy 



Duccio Fanelli 

Dipartimento di Energetica, Universita degli Studi di Firenze, 
via S. Maria 3, 50139 Florence, Italy and IN FN, Sezione di Firenze. 

, Francesca Di Patti 

■ Dipartimento di Fisica "Galileo Galilei", Universita degli Studi di Padova, via F. Marzolo 8, 35131 Padova, Italy 

rO ' A stochastic version of the Brusselator model is proposed and studied via the system size ex- 

pansion. The mean-field equations are derived and shown to yield to organized Turing patterns 
within a specific parameters region. When determining the Turing condition for instability, we pay 
particular attention to the role of cross diffusive terms, often neglected in the heuristic derivation 
of reaction diffusion schemes. Stochastic fluctuations are shown to give rise to spatially ordered 
solutions, sharing the same quantitative characteristic of the mean-field based Turing scenario, in 
term of excited wavelengths. Interestingly, the region of parameter yielding to the stochastic self- 
organization is wider than that determined via the conventional Turing approach, suggesting that 
the condition for spatial order to appear can be less stringent than customarily believed. 
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I. INTRODUCTION 
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Turing instability constitutes a universal paradigm for the spontaneous generation of spatially organized patterns 
[l[ . It formally applies to a wide category of phenomena, which can be modeled via the so called reaction diffusion 
schemes. These are mathematical models that describe the coupled evolution of spatially distributed species, as driven 
by microscopic reactions and freely diffusing in the embedding medium. Diffusion can potentially seed the instability 
by perturbing the mean-field homogeneous state, through an activator-inhibitor mechanism, and so yielding to the 
emergence of patched, spatially inhomogeneous, density distribution [2|. The realm of application of the Turing 
ideas encompasses different fields, ranging from chemistry to biology, from ecology to physics. The most astonishing 
■ examples, as already evidenced in Turing original paper, are perhaps encountered in the context of morphogenesis, 
the branch of embryology devoted to investigating the development of patterns and forms in biology Q . 

Beyond the qualitative agreement, one difficulty in establishing a quantitative link between theory and empirical 
' observations has to do with the strict conditions for which the organized Turing patterns are predicted to occur. In 
particular, and with reference to simple predator-prey competing populations, the relative degree of diffusivity of the 
interacting species has to be large according to the theory prescriptions, and at variance with the direct experimental 
evidence. Moreover, patterns formation appears to be rather robust in nature, as opposed to the Turing predictive 
scenario, where a fine tuning of the parameters is often necessary. 

In Q it was demonstrated that collective temporal oscillations can spontaneously emerge in a model of population 
dynamics, as due to a resonance mechanism that amplifies the unavoidable intrinsic noise, originating from the 
discreteness of the system. Later on an extension of the model was proposed Q so to explicitly account for the 
notion of space. It was in particulary shown that the number density of the interacting species oscillates both in 
time and space, a macroscopic effect resulting from the amplification of the stochastic fluctuations about the time 
independent solution of the deterministic equations. More recently, Butler and Goldenfeld [(| proved that persistent 
spatial patterns and temporal oscillations induced by demographic noise, can develop in a simple predator-prey model 
of plankton-herbivore dynamics. The model considered by the authors of [6( exhibits a Turing order in the mean 
field theory. The effect of intrinsic noise translates however into an enlargement of the parameter region yielding to 
the Turing mechanism, when compared to its homologous domain predicted within the conventional linear stability 
analysis. This is an individual based effect, which has to be accommodated for in any sensible model of natural 
phenomena, and which lacks in the Turing interpretative scenario, that formally applies to the idealized continuum 
limit. Interestingly, as reported in Q, the discreteness of the scrutinized medium can yield to robust spatio-temporal 



2 



structures, also when the system does not undergo Turing order in its mean-field, deterministic version. 

Starting from this setting, we present the results of our investigations carried out for a spatial version of the 
Brusselator model, which we shall be introducing in the forthcoming section. As opposed to the analysis in [f|, 
we will operate within the so called urn protocol, where individual elements belonging to the inspected species 
are assumed to populate an assigned container. The proposed formulation of the Brusselator differs from the one 
customarily reported in the literature. Additional terms are in fact obtained building on the underlying microscopic 
picture. These latter contributions are generally omitted, an assumption that we interpret as working in a diluted 
limit. Interestingly, the phase diagram of the homogeneous (aspatial) version of the model, is remarkably different 
from its diluted analogue, a fact we will substantiate in the following. 

Furthermore, when constructing the mean-field dynamics from the assigned microscopic rules, one recovers non- 
trivial cross diffusive terms, as previously remarked in These are derived within a self-consistent analysis and 
ultimately stands from assuming a finite carrying capacity in a given spatial patch. As such, they potentially bear 
an important physical meaning which deserves to be further elucidated. This scheme, alternative to conventional 
reaction-diffusion models, calls for an extension of the original Turing condition, that we will derive in the following. 
This is the second result of the paper. 

Finally, by making use of the van Kampen system size expansion, we are able to estimate the power spectrum of 
fluctuations analytically and so delineate the boundaries of the spatially ordered domain in the relevant parameters' 
space. As already remarked in Q, the role played by the intimate grainess can impact dramatically the Turing vision, 
returning a generalized scenario that holds promise to bridge the gap with observations. 

II. A SPATIAL VERSION OF THE BRUSSELATOR MODEL 

The Brusselator model represents a paradigmatic example of an autocatalytic chemical reaction. This scheme was 
originally devised in 1971 by Prigogine and Glandsdorff [8| and quickly gained its reputation as the prototype model 
for oscillating chemical reaction of the Belosouv-Zhabotinsky |9j, [l0( type. In the following we shall consider a slightly 
modified version of the original formulation, where the number of reactants X and Y is conserved and totals in N, 
including the empties, here called E Moreover we will consider a spatially extended system composed of Q, cells, 
each of size I, where reactions are supposed to occur Practically each cell hosts a replica of the Brusselator 

system: The molecules are however allowed to migrate between adjacent cells, which in turn implies an effective 
spatial coupling imputed to the microscopic molecular diffusion. A periodic geometry is also assumed so to restore 
the translational invariance. Although the calculation can be carried out in any space dimension D (see 0, 0|) we 
shall here mainly refer to the D = 1 case study, so to privilege the clarity of the message over technical complications. 
It should be however remarked that our conclusions are general and remain unchanged in extended spatial settings. 
Mathematically the model can be cast in the form: 

A + E t A A + X u 
Xi + B 4 Yi+B, 
2Xi + Yi 4 3A 4 , 
Xi — > Ei, 

where the index i runs from 1 to 57 and identifies the cell where the molecules are located. The autocatalytic species 
of interest to us are Xi and Yi. The elements A and B work as enzymatic activators, and keep constant in number. 
As such, they can be straightforwardly absorbed into the definition of the reaction rates. Let us label with rii (resp. 
nii and Oj) the number of element of type Xi (resp. Yi and Ei) populating cell i. Then, assuming N to identify the 
maximum number of available cases within each cell, one gets N = + m, + Oi a condition which can be exploited to 
reduce the actual number of dynamical variables to two. The dynamics of the model is ultimately related to studying 
the coupled interaction between the discrete species rii and mj. The migration between neighbors cells is specified 
through: 



Xi + Ej A Ei + Xj , 
Yi + Ej 4 Ei + Yj. 
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Assuming a perfect mixing in each individual cell i, the transition probabilities T(-|-) read 

N — n, ■ — m,i 



T(n,; + 1, rrii\rii, rrii) = a 



Nil 

T(m — l,mi + l\ni,mi) = b 



T(n,i + 1, nii - l\rii, rrii) = d , 

T(m - l,mi\rii,mi) = c j^-, 

where according to the standard convention, the rightmost input specifies the original state and the other entry stems 
for the final one. In addition, the migration mechanism between neighbors cell is controlled by: 

. m N — ri j — Too 
T{rii - l,7ij + ±\m,nj) = fi 



N NVz (2) 

m, Jy — rij — m 7 - 
T(to, - l,mj + l\mi,mj) = 8 — , 

where the positive constants /i and 8 quantify the diffusion ability of the two species and z stands for the number of 
first neighbors. To complete the notation setting, we introduce the fi-dimensional vectors n and m to identify the 
state of the system. Their i — th components respectively read n, and rrii. 

The aforementioned system is intrinsically stochastic. At time t there exists a finite probability to observe the 
system in the state characterized by n and m. Let us label P(n, m, t) such a probability. One can then write down 
the so called master equation, a differential equation which governs the dynamical evolution of the quantity P(n, m, t). 
The master equation for the case at hand takes the form 



d n 
-P(n,m,i) = £ 



(e x . - 1) T(rii + l,mi\ni,mi) + (e£. - 1) T(rii - I, rrii, K,mi) 



jei 



(3) 



+ ( e yi e y 3 _ ■*-) T{rrii — l,mj + l\mi,mj)] 



P(n,m,i), 



where use has been made of the definition of the step operators 

e%J{...,rii, ...,m) =/(..., n, ± 1, m), e^./(n, . . . , m u . . .) = /(n, . . . , m, ± 1, . . .), (4) 

and where the second sum in equation Q runs on first neighbors j G i. Equation Q is exact, the underlying 
dynamics being a Markov process. The master equation ([3]) contains information on both the ideal mean-field 
dynamics (formally recovered in the limit of diverging system size) and the finite N corrections. To bring into 
evidence those two components, one can proceed according to the prescriptions of van Kampen and write the 
normalized concentration relative to the interacting species as: 



where £j and to stand for the stochastic contribution. The ansatz ([S]) is motivated by the central limit theorem and 
holds provided the dynamics evolve far from the absorbing boundaries (extinction condition). Here 1/y/N plays the 
role of a small parameter and paves the way to a perturbative analysis of the master equation. At the leading order 
one recovers the mean-field equations, while the fluctuations are characterized as next to leading corrections. The 
following section is devoted to discussing the system dynamics, according to its mean-field approximation. 



III. THE MEAN-FIELD APPROXIMATION AND THE TURING INSTABILITY 



Performing the perturbative calculation one ends up with the following system of partial differential equations for 
the 4>i(t) and m cei l * : 
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d T 4>i = ~ Hi - d(j>i + a (1 - 4>i - ipi) + af%ipi + n [Afc + fcAipi - ipiAfc] , 
d T ipi = bfa - afriipi +5[Aipi + ipiAfa - fcAtpi] , 

where r = t/(NQ). In cq. ([5]) we have introduced the discrete Laplacian operator A acting on a generic function fi 
as 

A/i = !£(£-*)■ (7) 

jei 

The above sum runs on the first neighbors which we assume to total in z. The limit f2 — > oo corresponds to shrinking 
the lattice spacing I to zero and so obtaining the continuum mean-field description. In this limit the system © 
converges to: 

d T cj) = - b(j) - d<$> + a (1 - 4> - VO + c0V t + m[V 2 + 0V 2 i/> - VV 2 0], 

= 60 - aj> 2 ip + S[W 2 tp + ipV 2 (f> - <j>V 2 ip], ' S ! 

where the pop ulation fractions go over the population densities and the rescaling fi — > [il and S — > SI 2 have been 
performed jl4|. As also remarked in cross diffusive terms of the type <ft\7 2 ip — tp\7 2 (j) appear in the mean-field 
equations, as a relic of the microscopic rules of interaction. More precisely, they arise as a direct consequence of 
the finite carrying capacity hypothesis. This observation materializes in a crucial difference with respect to the 
heuristically proposed reaction diffusion schemes and points to the need for an extended Turing analysis, i.e. derive 
the generalized conditions yielding to organized spatial structures. Before addressing this specific issue, we start by 
discussing the homogeneous fixed point of ©. From hereon we will set a = d = 1, a choice already made in [l2| for 
the original Brusselator model and which will make possible to visualize our conclusion in the reference plan (6, c) for 
any fixed ratio of the diffusivity amount S and ji. 



A. The homogeneous fixed points 

Plugging a = d = 1 into equations ((SJ) and looking for homogeneous solutions implies 

(j> = (l -4> - ip) - (i + b)4> + c<j) 2 ijj, 



(9) 

I ip — b(f> — c<fi 2 ip. 

As an important remark we notice that in the diluted limit <f>, ip << 1, the system tends to the mean- field of 



the original Brusselator model, as e.g. reported in 12[. For a proper tuning of the chemical parameters (assigning 
significantly different strength) one can prove that, if initialized so to verify the diluted limit, the system stays 
diluted all along its subsequent evolution. As a corollary, the diluted solution is contained in the formulation of 
the Brusselator here considered, this latter being therefore regarded as a sound generalization of the former. Notice 
that when operating under diluted conditions, the cross diffusive terms appearing in the spatial equations flS]) can be 
neglected and the standard reaction-diffusion scheme recovered. Again, let us emphasize that it is the finite carrying 
capacity assumption that modifies the mean-field description. 
System Q admits three fixed points, namely: 




_ c-V-Sbc+e 2 (2 _ c+y/-8be+c 2 

- 4 c J 03 - 

J, _ 1 I V-8fcc+c 2 |.I _ 1 



x - .s/.r- ,■- ) ? _ 1 _ ■., -8bc+e 2 



The first of the above corresponds to the extinction of the species X: It exists for any choice of the freely changing 
parameters 6, c and corresponds to a stable attractor for the dynamics. This solution is not present in the original 
Brusselator mean-field equations and its origin can be traced back to the role played by the limited capacity of the 
container. 

The remaining two fixed points are respectively a saddle point and a (non trivial) stable attractor. They both 
manifest when the condition — 8& + c > is fulfilled. The saddle point partitions the available phase space into two 
domains, each defining the basin of attraction of the stable points. When c = 86 a saddle-node bifurcation occurs: the 
fixed points (02, ^2) and (03, V^) collide to eventually disappear. The finite carrying capacity that follows the "urn 
representation" of the Brusselator dynamics destroys the limit cycle solution which is instead found in its celebrated 
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classical analogue [12J: No periodic solutions are found in the mean-field approximation, unless the diluted limit is 
considered. The specific nature of fixed point (</> 3 , ip 3 ) changes as a function of the parameters (6, c), as illustrated 
in figure Q] In region I, it is a stable node, while in II it is a stable spiral. Region III identifies the parameters 
values for which it disappears. We are here concerned with the stability of the homogeneous fixed point (^3, ^3) to 
inhomogeneous perturbation: Can an instability develop and yield to organized Turing-like spatial patterns? 



c 




1 2 3 4 5 

FIG. 1: The nature of the fixed point (03, ^3) as a function of the parameters (b, c): In region I the fixed point is a stable node, 
while in II it is a stable spiral. The dashed line c = 86 delineates region III where the fixed point disappears, following a saddle 
node bifurcation. 



B. Extending the Turing instability mechanism: the effect of cross-diffusion 

To answer the previous question one has to perform a linear stability analysis around the selected homogeneous 
solution (03,V>3), hereafter termed (0, if>), accounting for the effect of non homogeneous perturbation. The formal 
development is closely inspired by the original Turing calculation, but now the effects of cross diffusive terms need 
to be accommodated for. The technical details of the derivation are enclosed in the annexed Appendix, where the 
most general problem is defined and inspected. We will here take advantage from the conclusion therein reached and 
present the results relative to the Brusselator model. In this specific case a = d = 1, the dispersion relation (|2"T)) reads 
in particular: 



A(fc 2 ) =A + k 2 B 



1 

16c 



C + Vk 2 + £k A 



(11) 



where 



A 



lib 



B 



1 



16 16 



36 ^fc{~8b + c)5 



8c 



4 



c)n 



4c 



C =256c 2 + 1286c 2 + 1446 2 c 2 - 32c 3 - 326c 3 + 2c 4 + y/c(-8b + c) (-32c 2 - 246c 2 + 2c 3 ) , 
V = - 64c 2 5 + 166c 2 6 + 8c 3 S + 128c 2 n + 966cV - 16cV + \/c(-86 + c)- 

• (-64c<5 - 80bc5 + 8c 2 6 + 128c^ + 326c^ - 16c 2 //) , 
£ = - 32bcS 2 + A0c 2 5 2 + I28bc6fi - 32c 2 S^i - 1286c/i 2 + 32cV + ^/c{-8b + c) (-24c5 2 - 32cdfi + 32c/x 2 



(12) 
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The perturbation gets amplified if A(fc 2 ) > over a finite interval of k values. The edges, fci and &2, of the interval 
are identified by imposing A(fc 2 ) = which returns the following results: 



A'i. 



where 



1 



2cS{b - 1) + 86 2 (/i - S) - bcfi + y/c(-8b + c)(26 - 2bS - bp,) ± (T + ^c(-8b + c) g) 1 



(13) 



T = 2(4c 2 <5 2 - 86c(2 + c)S 2 + 32b 4 (8 - p) 2 - 4b 3 c(85 2 - 25p + 3/t 2 ) + fe 2 c(4(12 + c)S 2 + cp 2 )), 

g = 2(-4cS 2 + 8bcS 2 + 8b 3 {25 2 -<*/*- p 2 ) + 6 2 (-4(4 + c)S 2 - 165/1 + cp 2 )). ^ 



By making use of conditions (|29|l as derived in the appendix, the predicted region for the Turing instability is 
traced in the (6, c) parameter space for an assigned ratio of mutual diffusivity, see figure O When compared to the 
conventional Turing analysis, based on the heuristically hypothesized reaction diffusion scheme (i.e. removing the 
cross terms in eq. l[5])). the region of inhomogeneous instability shrinks, the condition for spatial self-organization 
becoming even more peculiar than so far believed. 

Up to now we have carried out the analysis in the continuous mean-field scenario neglecting the role of finite size 
corrections. Stochastic fluctuations can however amplify due to an inherent resonant mechanism and consequently 
give rise to a large scale spatial order, similar to that predicted within the Turing interpretative picture Q- The 
following section is entirely devoted to clarifying this important point. 
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FIG. 2: Upper panel: The domain corresponding to the Turing instability are displayed in the (b, c) plan for 8 = 15 and /t = 1. 
The solid line (zone 2) delineates the region calculated when accounting for the cross-diffusive contribution, the dashed one 
(zone 1) stands for the conventional Turing analysis where diffusion is modeled via Laplacian operators (i.e. disregarding cross 
diffusion). Lower panel: Dispersion relation X(k) vs. k. Here c = 70 and b = 3, 5, 6, 7 (starting from the dotted bottom curve). 



7 



IV. THE ROLE OF STOCHASTIC NOISE AND THE POWER SPECTRUM OF FLUCTUATIONS 



At the next to leading order in the van Kampen perturbative development one ends up with a Fokker-Planck 
equation for the probability distribution of the fluctuations H(£,r),t) = P(n, m, t) where the vectors £ and rj have 
dimensions f2 and respective components ^ and r\i . The derivation is lengthy but straightforward and the reader can 
refer to e.g. [H, 0] for a detailed account on the mathematical technicalities. We will here report on the outcome of 
the calculation in terms of the predicted power spectra of the fluctuations close to equilibrium. In complete analogy 
with e.g. equation (36) and (37) of t 5| the latter can be written as: 

p KX (u) =< |4H| 2 > 



P fe ,rH =< \r] k {uj)\ 2 > 

where the u> and k stand respectively for Fourier temporal and spatial frequencies 15] . Exploiting the formal analogy 
between the governing Fokker-Planck equation and the equivalent Langevin equation one eventually obtains (see eqs. 
(38) and (39) in @): 



P k ,Y(uj) 



Ci+B fc , n ^ 2 
( w 2-na )2 + lV' 

C 2 + Bfc j22 W 2 



where: 



C k ,x(u) = B kA1 Ml 22 -2B kA2 M kA2 M ka2 +B ka2 Ml 12 , (15) 
h 11 - 2Bfc )12 Mfc i 2iM/i. jll + B fcjll M^ 21 . 



Cfc.yH = B fc!22 M 2 n -2B M2 M fei21 M M1 +B M1 Af 2 21 . (16) 




The 2x2 matrix M k of elements M k ,ij reads 

-2 - b + 2cM) + a ( 1 - if) A k -1 
M fe = I _ *\ V > . ( ■ \ " I , (17) 

b - 2c# + S^A k -c(j> 2 • ■ ' 

and the elements B^y respectively are: 

B fc ,n = 1 + b(f> - $ + c^tP - 2/i0Afc + 2/i0 2 A fe + 2/i#A fc , (18) 
B fe , 12 - -60-c0 2 ^, (19) 
B fe , 22 = b4> + cftiP - 254>A k + 2SHjA k + 28^ A k , (20) 

with the additional condition Bk,i% = Bfe. 2 i. Finally, i7 2 = detMfc, Tk — —tr'M.k and the Fourier transform of the 
Laplacian A k = 2[cos(kl) — 1]. In the continuum limit, the cell size goes to zero and scales as ~ k 2 

We are now in a position to represent the power spectrum of the fluctuations as predicted by the system size 
expansion. In figure [3] we display a selection of power spectra relative to one population and for different values of the 
parameter 6, at fixed c and diffusivity amount (see caption). The snapshots refer to a region of the parameter space for 
which we do not expect spatial order to appear, based on the Turing paradigm. However, and beyond the simplified 
mean-field viewpoint, a clear spatial peak is displayed, which gains in potency as the boundary of the Turing domain 
is being approached. Physically, it seems plausible to assume that the demographic noise can modify the dispersion 
relation. As a consequence, the curve A(fc), negative defined beyond the region of Turing instability, can locally cross 
the zero line, taking positive values in correspondence of specific k. These latter modes get therefore destabilized, 
yielding to quasi- Turing structures. Clearly, the k candidate to drive the stochastic instabilities are the ones close to 
kmax, the wave number that identifies the position of the maximum of the dispersion relation. If the above scenario 
is correct, k max should then be reasonably similar to the k value that locates the peak position in the power spectra 
P(k, -)(0). In figure (j4]) such comparison is drawn, for both species and over a window of b values, confirming the 
adequacy of the proposed interpretation. Based on the above, we can convincingly argue that the spatial modes here 
predicted represent an ideal continuation of the Turing structures beyond the region of parameters deputed to the 
mean-field instability and, for this reason, we suggest to refer to them as to stochastic Turing patterns. 
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The fact that spatial order appears for a wider range of the control parameter is further confirmed by visual 
inspection of figure [5j where the region of stochastic induced spatial organization, delineated from the power spectra 
calculated above, is depicted and compared to the corresponding mean-field prediction. Notice that the domains 
yielding to stochastic Turing patterns structures are different depending on the considered species. This observation 
marks yet another difference with respect to the standard mean-field theory. We can in fact expect to observe spatially 
ordered structures for just one of the two species, while the other is homogenously distributed. 

As an additional a point, we consider the case 5 = \i = 1. This choice cannot yield to the genuine Turing order, a 
pronounced difference in the relative diffusivity rates being an unavoidable pre-requisite. As opposed to the classical 
picture, stochastic Turing patterns can still emerge as demonstrated in figure [6] 

From the above expressions for the power spectra, one can also show that for fj, = 6 = (the aspatial limit) quasi- 
periodic time oscillations manifest in the concentrations. This is an effect of inner stochastic noise which restores the 
oscillations destroyed by the introduction od the carrying capacity. In general terms, self-oscillatory dynamics could 
be possibly understood, as resulting from the discrete nature of the medium, a scenario alternative to any ad hoc 
mean-field interpretation. 




12345 12345 



FIG. 3: Upper panels: The power spectra Pk,x(oo) for c = 70 and b = 3, 6 (going from left to right). Here 8 = 15 and (x = 1. 
The selected values of b and c fall outside the region of Turing instability, as follows the mean-field linear stability calculation 
(see figure[2j). Lower panels: The power spectra Pk,x(0) plotted vs. k so to enable one to appreciating the range of excited 
spatial wavelengths. This latter approximately matches that found in the realm of the Turing mean-field calculation (see figure 
[2}, implying similar characteristic of the spatially ordered structures. 



V. CONCLUSION 

Spatially organized patterns are reported to occur in a large gallery of widespread applications and are currently 
interpreted by resorting to the paradigmatic Turing picture. This vision, though successful, relies on a mean-field 
description of the relevant reaction diffusion schemes, often guessed on purely heuristic basis. An alternative scenario 
would require accounting for the intimate microscopic dynamics and so encapsulating the effect of the small scale 
grainess. This latter translates into stochastic fluctuations that, under specific conditions, can amplify and so give 
rise to organized spatio-temporal patterns. In this paper we have considered a microscopic model of Brusselator, 



9 




b 



4.0 4.5 5.0 5.5 6.0 6.5 7.0 4.0 4.5 5.0 5.5 6.0 6.5 7.0 



FIG. 4: Left panel: Symbols refer to the values in k that identify the maximum of the function Pk,x(0), here plotted as a 
function of b. The solid line stands for the position of the maximum of the mean-field dispersion relation X(k). Parameters are 
set as in figure [3] Right panel: as in the left panel, but now the symbols refer to species Y. 
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FIG. 5: The extended domain of Turing like instability predicted by the stochastic based analysis: the dashed line refers to 
species X, while the dot-dashed to species Y. The stochastic Turing region is confronted to the corresponding mean-field 
solution (solid line, also depicted in figure |3). 

P x (co=0, k) 
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1 2 3 4 5 



FIG. 6: Left panel: The power spectra Pfc,x(w) for a = d = 1, c = 70 and 6 = 8. Spatial and temporal peaks appear now to 
be decoupled. Right panel: The corresponding power spectra Pk,x(0) plotted vs. k. A clear peak is displayed. 
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and shown that Turing like patterns can indeed emerge beyond the parameter region predicted by the conventional 
Turing theory and due to the role of finite size corrections to the mean-field idealized dynamics. This result is 
obtained via a system size expansion which enables us to return closed analytical expressions for the power spectra of 
fluctuations. Our results agree with the conclusion reached in Q for another model and employing different analytical 
tools. Organized patterns can therefore occur more easily than expected, an observation that can potentially help 
reconciling theory and observations. In particular we find stochastic Turing patterns to emerge for S ~ fi, a condition 
for which Turing order is prevented to occur. 
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VI. APPENDIX 

Let us start by the generic set of partial differential equations: 

d t 4> = /(</>, t/>) + a*[v 2 + ^vV - ^ 2 <t>] , 
d t ifj = g (<f>, v>) + s [v 2 v> + V-vV - 0V 2 V>] , 

which explicitly allow for cross diffusive contributions. We choose to deal with zero flux boundary conditions. 

Assume that in absence of diffusion (homogeneous setting) the model tends towards a fixed point specified by (<f>, ?/>). 
In other words, (f>, ifi, do not depend over space variables. Allowing for diffusion, and imposing fi ^ 5, can make the 
system unstable to spatial perturbation. To clarify this point, let us start by considering the Jacobian matrix for 
H = 5 = Q: 



f<t> fip 
9<t> 9i> 



(22) 



The fixed point (</>, ip) is linearly stable if J has positive determinant and negative trace: 

det J = fyg^ - g$f$ > 0, trJ = fa + g^ < 0. (23) 

Assuming (|2"31 to hold one can proceed with a linearization of (|2"Tj) (with fi, 6 ^ 0) around (<f>, ip) and so looking for 
the sought condition of instability. Define x = then (|2Tj) become 

7 2„ ™ _ ( v C 1 - 



a «x^x + DVV D^-^-.J. ,24) 

Define the eigenfunctions of the Laplacian operator as: 

(V 2 + fc 2 ) W fe (r) = 0, re If, 
and write the solution to eq. (|24[) in the form: 

x(t,r) =Y / e Xt a k W k (r). (25) 

k 

Assume fj, = 1, or equivalently label with 5 the ratio of the two diffusivities (after proper rescaling of the rate 
coefficients). Substituting ansatz (|2"5)) into eq. §Zty yields: 

e At [j-fc 2 D- Al]W fe =0, 

which implies that the system admits a solution iff the matrix J — k 2 D — Al is singular, i.e.: 
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det(J — fc 2 D — Al) = 0. (26) 

Label the solutions of (|26p as Ai(fc 2 ) and A2(fc 2 ): They can be interpreted as dispersion relation, specifying the 
time scale of departure (or convergence) of the fc-th mode towards the deputed fixed point. If at least one of the 
two solutions displays a positive real part, the mode is unstable, and drives the system dynamics towards a non- 
homogeneous configuration in response to the initial perturbation. Manipulating the determinant, (|26[) takes the 
form: 



X 2 + X q(k 2 ) + h(k 2 ) = 0, (27) 

where 

q(k 2 ) =k 2 + k 2 5(l -$-$)- fy- 

h(k 2 ) = k 4 5 -k 4 4>5- k 4 4> S - fc 2 <S/ + k 2 4> &U + fe2 ^ S U + fc 2 050" 

U 94> - k2 9i> + k2 ^P 9i> + f<p9i>- 

Xi(k 2 ) and X2(k 2 ) are obtained as roots of (f2"7| . The actual dispersion relation, X(k 2 ) in the main body of the paper, 
is the one that displays the largest real part. We notice that (i) q(k 2 ) is always positive as 6 > by definition; (ii) 
(1 — <j> — if)) > as and if; are concentrations; (iii) —f$ — g^ > due to (|23|) . Left-hand side of equation ((|27p) 
is hence a parabola with the concavity pointing upward and the minimum positioned in the half-plane of negative 
abscissa. In conclusion a necessary and sufficient condition for the existence of real roots is h < where: 

[a =5(l-4>-4>), 

h{k 2 ) = ~ak 4 + bk 2 + d, lb =4/ r ^ + ^ + ^) + ^+ ft ), (28) 

Ic =detJ, 

with a > 0, c > 0. The condition for ft, < corresponds to b < and b — 4ac > 0. Summing up the condition for the 
generalized Turing instability reads: 

(SU + gi,)-(j) (6 fa +g4,)-i> (SU + g^) > 0, 

2 (29 J 

[tfU + ffv) - <t> + 9<f>)-^ (SU + 9i>)] > 4 (5(1 - - V) det J, 

together with (|23|) . 
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